#setwd('/afs/inf.ed.ac.uk/user/s17/s1725186/Documents/PhD-Models/FirstPUModel/RMarkdowns')

library(tidyverse) ; library(reshape2) ; library(glue) ; library(plotly) ; library(dendextend)
library(RColorBrewer) ; library(viridis) ; require(gridExtra) ; library(GGally)
library(expss)
library(polycor)
library(foreach) ; library(doParallel)
library(knitr)
library(biomaRt)
library(anRichment) ; library(BrainDiseaseCollection)
suppressWarnings(suppressMessages(library(WGCNA)))

SFARI_colour_hue = function(r) {
  pal = c('#FF7631','#FFB100','#E8E328','#8CC83F','#62CCA6','#59B9C9','#b3b3b3','#808080','gray','#d9d9d9')[r]
}

Load preprocessed dataset (preprocessing code in 19_10_14_data_preprocessing.Rmd) and clustering (pipeline in 19_10_21_WGCNA.Rmd)

# Gandal dataset
load('./../Data/preprocessed_data.RData')
datExpr = datExpr %>% data.frame
DE_info = DE_info %>% data.frame


# GO Neuronal annotations: regex 'neuron' in GO functional annotations and label the genes that make a match as neuronal
GO_annotations = read.csv('./../Data/genes_GO_annotations.csv')
GO_neuronal = GO_annotations %>% filter(grepl('neuron', go_term)) %>% 
              mutate('ID'=as.character(ensembl_gene_id)) %>% 
              dplyr::select(-ensembl_gene_id) %>% distinct(ID) %>%
              mutate('Neuronal'=1)


# SFARI Genes
SFARI_genes = read_csv('./../../../SFARI/Data/SFARI_genes_08-29-2019_w_ensembl_IDs.csv')
SFARI_genes = SFARI_genes[!duplicated(SFARI_genes$ID) & !is.na(SFARI_genes$ID),]


# Clusterings
clusterings = read_csv('./../Data/clusters.csv')


# Update DE_info with SFARI and Neuronal information
genes_info = DE_info %>% mutate('ID'=rownames(.)) %>% left_join(SFARI_genes, by='ID') %>% 
  mutate(`gene-score`=ifelse(is.na(`gene-score`), 'None', `gene-score`)) %>%
  left_join(GO_neuronal, by='ID') %>% left_join(clusterings, by='ID') %>%
  mutate(Neuronal=ifelse(is.na(Neuronal), 0, Neuronal)) %>%
  mutate(gene.score=ifelse(`gene-score`=='None' & Neuronal==1, 'Neuronal', `gene-score`), 
         significant=padj<0.05 & !is.na(padj))

# Add gene symbol
getinfo = c('ensembl_gene_id','external_gene_id')
mart = useMart(biomart='ENSEMBL_MART_ENSEMBL', dataset='hsapiens_gene_ensembl',
               host='feb2014.archive.ensembl.org') ## Gencode v19
gene_names = getBM(attributes=getinfo, filters=c('ensembl_gene_id'), values=genes_info$ID, mart=mart)

genes_info = genes_info %>% left_join(gene_names, by=c('ID'='ensembl_gene_id'))


clustering_selected = 'DynamicHybrid'
genes_info$Module = genes_info[,clustering_selected]

dataset = read.csv(paste0('./../Data/dataset_', clustering_selected, '.csv'))
dataset$Module = dataset[,clustering_selected]


rm(DE_info, GO_annotations, clusterings, getinfo, mart, dds)


Relation to external clinical traits


Quantifying module-trait associations


Using the hetcor function, that calculates Pearson, polyserial or polychoric correlations depending on the type of variables involved.

datTraits = datMeta %>% dplyr::select(Diagnosis, Region, Sex, Age, PMI, RNAExtractionBatch) %>%
            dplyr::rename('ExtractionBatch' = RNAExtractionBatch)

# Recalculate MEs with color labels
ME_object = datExpr %>% t %>% moduleEigengenes(colors = genes_info$Module)
MEs = orderMEs(ME_object$eigengenes)

# Calculate correlation between eigengenes and the traits and their p-values
moduleTraitCor = MEs %>% apply(2, function(x) hetcor(x, datTraits)$correlations[1,-1]) %>% t
rownames(moduleTraitCor) = colnames(MEs)
colnames(moduleTraitCor) = colnames(datTraits)
moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nrow(datExpr))

# Create text matrix for the Heatmap
textMatrix = paste0(signif(moduleTraitCor, 2), ' (', signif(moduleTraitPvalue, 1), ')')
dim(textMatrix) = dim(moduleTraitCor)

# In case there are any NAs
if(sum(!complete.cases(moduleTraitCor))>0){
  print(paste0(sum(is.na(moduleTraitCor)),' correlation(s) could not be calculated')) 
}
## [1] "1 correlation(s) could not be calculated"
rm(ME_object)

Note: The correlation between Module #F464E5 and Diagonsis is the one that cannot be calculated, weirdly enough, the thing that causes the error is that the initial correlation is too high, so it would be a very bad thing to lose this module because of this numerical error. I’m going to fill in its value using the polyserial function, which doesn’t give exactly the same results as the hetcor() function, but it’s quite similar.

# Calculate the correlation tha failed with hetcor()
moduleTraitCor['ME#F464E5','Diagnosis'] = polyserial(MEs[,'ME#F464E5'], datTraits$Diagnosis)
## Warning in polyserial(MEs[, "ME#F464E5"], datTraits$Diagnosis): initial
## correlation inadmissible, 1.01199040044846, set to 0.9999

I’m going to select all the modules that have an absolute correlation higher than 0.9 with Diagnosis to study them

# Sort moduleTraitCor by Diagnosis
moduleTraitCor = moduleTraitCor[order(moduleTraitCor[,1], decreasing=TRUE),]
moduleTraitPvalue = moduleTraitPvalue[order(moduleTraitCor[,1], decreasing=TRUE),]

# Create text matrix for the Heatmap
textMatrix = paste0(signif(moduleTraitCor, 2), ' (', signif(moduleTraitPvalue, 1), ')')
dim(textMatrix) = dim(moduleTraitCor)


labeledHeatmap(Matrix = moduleTraitCor, xLabels = names(datTraits), yLabels =  gsub('ME','',rownames(moduleTraitCor)), 
               yColorWidth=0, colors = brewer.pal(11,'PiYG'), bg.lab.y = gsub('ME','',rownames(moduleTraitCor)),
               textMatrix = textMatrix, setStdMargins = FALSE, cex.text = 0.8, cex.lab.y = 0.75, zlim = c(-1,1),
               main = paste('Module-Trait relationships'))

diagnosis_cor = data.frame('Module' = gsub('ME','',rownames(moduleTraitCor)),
                           'MTcor' = moduleTraitCor[,'Diagnosis'],
                           'MTpval' = moduleTraitPvalue[,'Diagnosis'])

genes_info = genes_info %>% left_join(diagnosis_cor, by='Module')

rm(moduleTraitPvalue, datTraits, textMatrix, diagnosis_cor)


Studying the modules with the highest absolute correlation to Diagnosis


The modules consist mainly of points with very high (absolute) values in PC2 (which we know is related to lfc), so this result is consistent with the high correlation between Module and Diagnosis, although some of the points with the highest PC2 values do not belong to these top modules

top_modules = gsub('ME','',rownames(moduleTraitCor)[abs(moduleTraitCor[,'Diagnosis'])>0.9])

cat(paste0('Top modules selected: ', paste(top_modules, collapse=', '),'\n'))
## Top modules selected: #F464E5, #00BECA
pca = datExpr %>% prcomp

plot_data = data.frame('ID'=rownames(datExpr), 'PC1' = pca$x[,1], 'PC2' = pca$x[,2]) %>%
            left_join(dataset, by='ID') %>% left_join(genes_info %>% dplyr::select(ID, external_gene_id), by='ID') %>%
            dplyr::select(ID, external_gene_id, PC1, PC2, Module, gene.score) %>%
            mutate(ImportantModules = ifelse(Module %in% top_modules, as.character(Module), 'Others')) %>%
            mutate(color = ifelse(ImportantModules=='Others','gray',ImportantModules),
                   alpha = ifelse(ImportantModules=='Others', 0.2, 0.4),
                   gene_id = paste0(ID, ' (', external_gene_id, ')'))

table(plot_data$ImportantModules)
## 
## #00BECA #F464E5  Others 
##    1638     626   13884
ggplotly(plot_data %>% ggplot(aes(PC1, PC2, color=ImportantModules)) + 
         geom_point(alpha=plot_data$alpha, color=plot_data$color, aes(ID=gene_id)) + theme_minimal() + 
           ggtitle('Modules with strongest relation to Diagnosis'))
rm(pca)




Module Membership vs Gene Significance


create_plot = function(module){
  
  plot_data = dataset %>% dplyr::select(ID, paste0('MM.',gsub('#','',module)), GS, gene.score) %>% filter(dataset$Module==module)
  colnames(plot_data)[2] = 'Module'
  
  SFARI_colors = as.numeric(names(table(as.character(plot_data$gene.score)[plot_data$gene.score!='None'])))
  
  p = ggplotly(plot_data %>% ggplot(aes(Module, GS, color=gene.score)) + geom_point(alpha=0.5, aes(ID=ID)) + ylab('Gene Significance') +
               scale_color_manual(values=SFARI_colour_hue(r=c(SFARI_colors,8))) + theme_minimal() + xlab('Module Membership') +
               ggtitle(paste0('Module ', module,'  (MTcor = ', round(moduleTraitCor[paste0('ME',module),1],2),')')))
  
  return(p)
}

create_plot(top_modules[1])
create_plot(top_modules[2])
rm(create_plot)




SFARI Genes


List of SFARI Genes in top modules ordered by SFARI score and Gene Significance

table_data = dataset %>% left_join(genes_info %>% dplyr::select(ID, external_gene_id), by='ID') %>%
             dplyr::select(ID, external_gene_id, GS, gene.score, Module) %>% arrange(gene.score, desc(abs(GS))) %>%
             dplyr::rename('Ensembl ID'=ID, 'Gene Symbol'=external_gene_id, 
                    'SFARI score'=gene.score, 'Gene Significance'=GS)


kable(table_data %>% filter(Module == top_modules[1] & `SFARI score` != 'None') %>% dplyr::select(-Module),
      caption=paste0('SFARI Genes for Module ', top_modules[1]))
SFARI Genes for Module #F464E5
Ensembl ID Gene Symbol Gene Significance SFARI score
ENSG00000038382 TRIO 0.6264110 2
ENSG00000168769 TET2 0.8178650 3
ENSG00000083168 KAT6A 0.5998865 3
ENSG00000113742 CPEB4 0.3701387 3
ENSG00000145020 AMT 0.3424890 3
ENSG00000132510 KDM6B 0.3330706 3
ENSG00000128573 FOXP2 0.2489053 3
ENSG00000112902 SEMA5A 0.2191666 3
ENSG00000165995 CACNB2 0.1701529 3
ENSG00000137801 THBS1 0.8076180 4
ENSG00000106366 SERPINE1 0.7907461 4
ENSG00000011422 PLAUR 0.7753875 4
ENSG00000172554 SNTG2 0.7354320 4
ENSG00000181481 RNF135 0.6768553 4
ENSG00000152208 GRID2 0.6205857 4
ENSG00000227184 EPPK1 0.5946476 4
ENSG00000183098 GPC6 0.5839280 4
ENSG00000188785 ZNF548 0.5801457 4
ENSG00000198589 LRBA 0.5579943 4
ENSG00000182372 CLN8 0.3972145 4
ENSG00000071246 VASH1 0.3828434 4
ENSG00000196277 GRM7 0.3773957 4
ENSG00000239389 PCDHA13 0.3419447 4
ENSG00000007314 SCN4A 0.3327403 4
ENSG00000182472 CAPN12 0.3178125 4
ENSG00000224389 C4B 0.3077509 4
ENSG00000215045 GRID2IP 0.2888427 4
ENSG00000130508 PXDN 0.1235447 4
ENSG00000187957 DNER 0.0565791 4
ENSG00000100024 UPB1 0.3133084 5
ENSG00000026508 CD44 0.2637995 5
ENSG00000136244 IL6 0.1960508 5
ENSG00000184588 PDE4B 0.1947772 5
ENSG00000141526 SLC16A3 0.8254828 6
ENSG00000072364 AFF4 0.7711840 6
ENSG00000171791 BCL2 0.3611981 6
kable(table_data %>% filter(Module == top_modules[2] & `SFARI score` != 'None') %>% dplyr::select(-Module),
      caption=paste0('SFARI Genes for Module ', top_modules[2]))
SFARI Genes for Module #00BECA
Ensembl ID Gene Symbol Gene Significance SFARI score
ENSG00000136535 TBR1 -0.6544304 1
ENSG00000136531 SCN2A -0.5379355 1
ENSG00000174469 CNTNAP2 -0.7002325 2
ENSG00000061676 NCKAP1 -0.6918489 2
ENSG00000139613 SMARCC2 -0.4332840 2
ENSG00000196557 CACNA1H -0.4139346 2
ENSG00000119866 BCL11A -0.3842245 2
ENSG00000157445 CACNA2D3 -0.3391182 2
ENSG00000144619 CNTN4 -0.3111072 2
ENSG00000074590 NUAK1 -0.8234159 3
ENSG00000144285 SCN1A -0.8004590 3
ENSG00000196876 SCN8A -0.7890809 3
ENSG00000170579 DLGAP1 -0.7792697 3
ENSG00000005955 GGNBP2 -0.7656571 3
ENSG00000132294 EFR3A -0.7504244 3
ENSG00000078328 RBFOX1 -0.7315557 3
ENSG00000003147 ICA1 -0.7117475 3
ENSG00000197535 MYO5A -0.7116705 3
ENSG00000175497 DPP10 -0.7061637 3
ENSG00000182621 PLCB1 -0.6921795 3
ENSG00000171759 PAH -0.6897098 3
ENSG00000021645 NRXN3 -0.6713796 3
ENSG00000166501 PRKCB -0.6155044 3
ENSG00000183454 GRIN2A -0.5294268 3
ENSG00000146830 GIGYF1 -0.5256914 3
ENSG00000185345 PARK2 -0.4980108 3
ENSG00000164506 STXBP5 -0.4858781 3
ENSG00000132024 CC2D1A -0.4128981 3
ENSG00000050030 KIAA2022 -0.4002974 3
ENSG00000168116 KIAA1586 -0.3735057 3
ENSG00000133026 MYH10 -0.3267293 3
ENSG00000101489 CELF4 -0.3181586 3
ENSG00000185008 ROBO2 -0.3077778 3
ENSG00000182256 GABRG3 -0.2873575 3
ENSG00000139174 PRICKLE1 -0.2147716 3
ENSG00000138411 HECW2 -0.1373931 3
ENSG00000140945 CDH13 0.0845059 3
ENSG00000166147 FBN1 -0.0824793 3
ENSG00000171723 GPHN -0.0749364 3
ENSG00000149972 CNTN5 -0.0485789 3
ENSG00000144406 UNC80 -0.8620647 4
ENSG00000135631 RAB11FIP5 -0.8335775 4
ENSG00000132639 SNAP25 -0.8106671 4
ENSG00000163618 CADPS -0.8086094 4
ENSG00000159082 SYNJ1 -0.7858980 4
ENSG00000104093 DMXL2 -0.7629031 4
ENSG00000172260 NEGR1 -0.7083293 4
ENSG00000163399 ATP1A1 -0.6896505 4
ENSG00000115840 SLC25A12 -0.6828765 4
ENSG00000196090 PTPRT -0.6742079 4
ENSG00000130477 UNC13A -0.6667346 4
ENSG00000154263 ABCA10 -0.6407178 4
ENSG00000153575 TUBGCP5 -0.6283954 4
ENSG00000136928 GABBR2 -0.6252262 4
ENSG00000144290 SLC4A10 -0.5769256 4
ENSG00000144331 ZNF385B -0.5383310 4
ENSG00000186094 AGBL4 -0.5377404 4
ENSG00000107105 ELAVL2 -0.5374477 4
ENSG00000196433 ASMT -0.5323119 4
ENSG00000155961 RAB39B -0.5256637 4
ENSG00000116141 MARK1 -0.5062338 4
ENSG00000155886 SLC24A2 -0.4946020 4
ENSG00000109158 GABRA4 -0.4943945 4
ENSG00000198010 DLGAP2 -0.4942182 4
ENSG00000243232 PCDHAC2 -0.4777463 4
ENSG00000128594 LRRC4 -0.4719222 4
ENSG00000175161 CADM2 -0.4567470 4
ENSG00000182901 RGS7 -0.4466930 4
ENSG00000133019 CHRM3 -0.4273514 4
ENSG00000172915 NBEA -0.4149285 4
ENSG00000156113 KCNMA1 -0.4061200 4
ENSG00000081189 MEF2C -0.4000890 4
ENSG00000114374 USP9Y -0.3520447 4
ENSG00000174442 ZWILCH -0.3276604 4
ENSG00000136653 RASSF5 -0.3267914 4
ENSG00000151790 TDO2 -0.3195670 4
ENSG00000152495 CAMK4 -0.3178490 4
ENSG00000060140 STYK1 -0.3136339 4
ENSG00000064687 ABCA7 -0.3112037 4
ENSG00000165355 FBXO33 -0.3050639 4
ENSG00000135127 CCDC64 -0.3047331 4
ENSG00000170289 CNGB3 -0.2731444 4
ENSG00000139915 MDGA2 -0.2430880 4
ENSG00000152402 GUCY1A2 -0.2327724 4
ENSG00000111615 KRR1 -0.2248392 4
ENSG00000100359 SGSM3 -0.2075939 4
ENSG00000159640 ACE -0.1924407 4
ENSG00000004468 CD38 -0.1674841 4
ENSG00000128513 POT1 -0.1434731 4
ENSG00000244588 RAD21L1 -0.1256607 4
ENSG00000165379 LRFN5 -0.1218335 4
ENSG00000100038 TOP3B -0.0993690 4
ENSG00000152214 RIT2 0.0792829 4
ENSG00000165246 NLGN4Y -0.0675768 4
ENSG00000181036 FCRL6 0.0604928 4
ENSG00000007168 PAFAH1B1 -0.8564947 5
ENSG00000087470 DNM1L -0.8491182 5
ENSG00000104725 NEFL -0.7912259 5
ENSG00000022355 GABRA1 -0.7777200 5
ENSG00000118596 SLC16A7 -0.7640017 5
ENSG00000171735 CAMTA1 -0.7302306 5
ENSG00000146469 VIP -0.6852620 5
ENSG00000152969 JAKMIP1 -0.6814938 5
ENSG00000106123 EPHB6 -0.6655863 5
ENSG00000173726 TOMM20 -0.6303781 5
ENSG00000100362 PVALB -0.4940183 5
ENSG00000007171 NOS2 -0.4656878 5
ENSG00000102468 HTR2A -0.4562790 5
ENSG00000145979 TBC1D7 -0.4227883 5
ENSG00000130707 ASS1 -0.3913746 5
ENSG00000158258 CLSTN2 -0.3471046 5
ENSG00000196730 DAPK1 -0.2855914 5
ENSG00000166603 MC4R -0.2783239 5
ENSG00000185666 SYN3 -0.2690820 5
ENSG00000157168 NRG1 -0.1662121 5
ENSG00000166866 MYO1A -0.1473913 5
ENSG00000106113 CRHR2 -0.0441381 5
ENSG00000077279 DCX -0.5750475 6
ENSG00000148680 HTR7 -0.4261417 6
ENSG00000125780 TGM3 -0.3834141 6

Modules with the strongest module-diagnosis correlation should have the highest percentage of SFARI Genes, but this doesn’t seem to be the case (even the opposite may be true)

plot_data = dataset %>% mutate('hasSFARIscore' = gene.score!='None') %>% 
            group_by(Module, MTcor, hasSFARIscore) %>% summarise(p=n()) %>% 
            left_join(dataset %>% group_by(Module) %>% summarise(n=n()), by='Module') %>% 
            mutate(p=round(p/n*100,2)) 

for(i in 1:nrow(plot_data)){
  this_row = plot_data[i,]
  if(this_row$hasSFARIscore==FALSE & this_row$p==100){
    new_row = this_row
    new_row$hasSFARIscore = TRUE
    new_row$p = 0
    plot_data = plot_data %>% rbind(new_row)
  }
}

plot_data = plot_data %>% filter(hasSFARIscore==TRUE)

ggplotly(plot_data %>% ggplot(aes(MTcor, p, size=n)) + geom_smooth(color='gray', se=FALSE) +
         geom_point(color=plot_data$Module, alpha=0.5, aes(id=Module)) + geom_hline(yintercept=mean(plot_data$p), color='gray') +
         xlab('Module-Diagnosis correlation') + ylab('% of SFARI genes') +
         theme_minimal() + theme(legend.position = 'none'))
rm(i, this_row, new_row, plot_data)




Module Eigengenes


Since these modules have the strongest relation to autism, this pattern should be reflected in their model eigengenes, having two different behaviours for the samples corresponding to autism and the ones corresponding to control.

In both cases, the Eigengenes separate the behaviour between autism and control samples very clearly!

plot_EGs = function(module){

  plot_data = data.frame('ID' = rownames(MEs), 'MEs' = MEs[,paste0('ME',module)], 'Diagnosis' = datMeta$Diagnosis)

  p = plot_data %>% ggplot(aes(Diagnosis, MEs, fill=Diagnosis)) + geom_boxplot() + theme_minimal() + theme(legend.position='none') +
                    ggtitle(paste0('Module ', module, '  (MTcor=',round(moduleTraitCor[paste0('ME',module),1],2),')'))
  return(p)
}

p1 = plot_EGs(top_modules[1])
p2 = plot_EGs(top_modules[2])

grid.arrange(p1, p2, nrow=1)

rm(plot_EGs, p1, p2)




Identifying important genes


Selecting the modules with the highest correlation to Diagnosis, and, from them, the genes with the highest module membership-(absolute) gene significance

*Ordered by \(\frac{MM+|GS|}{2}\)

There aren’t that many SFARI genes in the top genes of the modules and not a single one belonging to scores 1 and 2

create_table = function(module){
  top_genes = dataset %>% left_join(genes_info %>% dplyr::select(ID, external_gene_id), by='ID') %>% 
              dplyr::select(ID, external_gene_id, paste0('MM.',gsub('#','',module)), GS, gene.score) %>%
              filter(dataset$Module==module) %>% dplyr::rename('MM' = paste0('MM.',gsub('#','',module))) %>% 
              mutate(importance = (MM+abs(GS))/2) %>% arrange(by=-importance) %>% top_n(20)
  return(top_genes)
}

top_genes_1 = create_table(top_modules[1])
kable(top_genes_1, caption=paste0('Top 10 genes for module ', top_modules[1], '  (MTcor = ',
                                  round(moduleTraitCor[paste0('ME',top_modules[1]),1],2),')'))
Top 10 genes for module #F464E5 (MTcor = 1)
ID external_gene_id MM GS gene.score importance
ENSG00000143384 MCL1 0.8476929 0.9075109 None 0.8776019
ENSG00000161638 ITGA5 0.8510043 0.8600246 None 0.8555144
ENSG00000089159 PXN 0.8085675 0.8694821 None 0.8390248
ENSG00000196935 SRGAP1 0.7501468 0.9223277 None 0.8362373
ENSG00000003402 CFLAR 0.7941123 0.8169199 None 0.8055161
ENSG00000148841 ITPRIP 0.7881631 0.8162584 None 0.8022107
ENSG00000124782 RREB1 0.7589325 0.8451372 None 0.8020348
ENSG00000106366 SERPINE1 0.8065142 0.7907461 4 0.7986301
ENSG00000150457 LATS2 0.7443403 0.8501398 None 0.7972400
ENSG00000138119 MYOF 0.7514562 0.8343142 None 0.7928852
ENSG00000073792 IGF2BP2 0.7392534 0.8389676 None 0.7891105
ENSG00000138166 DUSP5 0.8112245 0.7636771 None 0.7874508
ENSG00000072364 AFF4 0.7963444 0.7711840 6 0.7837642
ENSG00000120278 PLEKHG1 0.7556018 0.8087123 None 0.7821570
ENSG00000162745 OLFML2B 0.7023713 0.8526646 None 0.7775179
ENSG00000154640 BTG3 0.7244112 0.8301676 None 0.7772894
ENSG00000133639 BTG1 0.7423866 0.8079258 None 0.7751562
ENSG00000168769 TET2 0.7312714 0.8178650 3 0.7745682
ENSG00000120690 ELF1 0.7613055 0.7876526 None 0.7744790
ENSG00000158615 PPP1R15B 0.7289451 0.8197298 None 0.7743375
top_genes_2 = create_table(top_modules[2])
kable(top_genes_2, caption=paste0('Top 10 genes for module ', top_modules[2], '  (MTcor = ',
                                  round(moduleTraitCor[paste0('ME',top_modules[2]),1],2),')'))
Top 10 genes for module #00BECA (MTcor = -0.92)
ID external_gene_id MM GS gene.score importance
ENSG00000050748 MAPK9 0.9055266 -0.9401065 None 0.9228165
ENSG00000108395 TRIM37 0.9056872 -0.9297206 None 0.9177039
ENSG00000138078 PREPL 0.8926949 -0.9043456 None 0.8985203
ENSG00000177432 NAP1L5 0.8738603 -0.9048896 None 0.8893750
ENSG00000163577 EIF5A2 0.8622485 -0.9110526 None 0.8866505
ENSG00000128881 TTBK2 0.8938083 -0.8754731 None 0.8846407
ENSG00000176490 DIRAS1 0.8610032 -0.8968398 None 0.8789215
ENSG00000171132 PRKCE 0.8477771 -0.9090963 None 0.8784367
ENSG00000155097 ATP6V1C1 0.8504461 -0.8966833 None 0.8735647
ENSG00000196876 SCN8A 0.9470763 -0.7890809 3 0.8680786
ENSG00000114573 ATP6V1A 0.8139207 -0.9175213 None 0.8657210
ENSG00000111674 ENO2 0.8761932 -0.8485455 None 0.8623694
ENSG00000125814 NAPB 0.8922318 -0.8313011 None 0.8617664
ENSG00000184368 MAP7D2 0.8332534 -0.8807598 None 0.8570066
ENSG00000131437 KIF3A 0.8671455 -0.8464587 None 0.8568021
ENSG00000144285 SCN1A 0.9127631 -0.8004590 3 0.8566110
ENSG00000162694 EXTL2 0.8023344 -0.8984415 None 0.8503880
ENSG00000172348 RCAN2 0.7828069 -0.9152892 None 0.8490480
ENSG00000130540 SULT4A1 0.8830724 -0.8116827 None 0.8473775
ENSG00000132639 SNAP25 0.8784733 -0.8106671 4 0.8445702
rm(create_table)
pca = datExpr %>% prcomp

plot_data = data.frame('ID'=rownames(datExpr), 'PC1' = pca$x[,1], 'PC2' = pca$x[,2]) %>%
            left_join(dataset, by='ID') %>% dplyr::select(ID, PC1, PC2, Module, gene.score) %>%
            mutate(color = ifelse(Module %in% top_modules, as.character(Module), 'gray')) %>%
            mutate(alpha = ifelse(color %in% top_modules & 
                                  ID %in% c(as.character(top_genes_1$ID), 
                                            as.character(top_genes_2$ID)), 1, 0.1))

plot_data %>% ggplot(aes(PC1, PC2)) + geom_point(alpha=plot_data$alpha, color=plot_data$color) + 
              theme_minimal() + ggtitle('Important genes identified through WGCNA')




Enrichment Analysis


Using the package anRichment

# Prepare dataset

# Create dataset with top modules membership and removing the genes without an assigned module
EA_dataset = data.frame('ensembl_gene_id' = genes_info$ID,
                        module = ifelse(genes_info$Module %in% top_modules, genes_info$Module, 'other')) %>%
             filter(genes_info$Module!='gray')

# Assign Entrez Gene Id to each gene
getinfo = c('ensembl_gene_id','entrezgene')
mart = useMart(biomart='ENSEMBL_MART_ENSEMBL', dataset='hsapiens_gene_ensembl', host='feb2014.archive.ensembl.org')
biomart_output = getBM(attributes=getinfo, filters=c('ensembl_gene_id'), values=EA_dataset$ensembl_gene_id, mart=mart)
## Cache found
EA_dataset = EA_dataset %>% left_join(biomart_output, by='ensembl_gene_id')

for(tm in top_modules){
  cat(paste0('\n',sum(EA_dataset$module==tm & is.na(EA_dataset$entrezgene)), ' genes from top module ',
             tm, ' don\'t have an Entrez Gene ID')) 
}
## 
## 8 genes from top module #F464E5 don't have an Entrez Gene ID
## 29 genes from top module #00BECA don't have an Entrez Gene ID
rm(getinfo, mart, biomart_output, tm)
# Manual: https://horvath.genetics.ucla.edu/html/CoexpressionNetwork/GeneAnnotation/Tutorials/anRichment-Tutorial1.pdf

collectGarbage()

# EA_dataset = rbind(EA_dataset[EA_dataset$module!='other',], EA_dataset[EA_dataset$module=='other',][sample(sum(EA_dataset$module=='other'), 1000),])

# Prepare datasets
GO_col = buildGOcollection(organism = 'human', verbose = 0)
## Loading required package: org.Hs.eg.db
## 
## 'select()' returned 1:many mapping between keys and columns
## 'select()' returned 1:1 mapping between keys and columns
internal_col = internalCollection(organism = 'human')
MillerAIBS_col = MillerAIBSCollection(organism = 'human')
BrainDisease_col = BrainDiseaseCollection(organism = 'human')
combined_col = mergeCollections(GO_col, internal_col, MillerAIBS_col, BrainDisease_col)

# Print collections used
cat('Using collections: ')
## Using collections:
knownGroups(combined_col, sortBy = 'size')
##  [1] "GO"                                         
##  [2] "GO.BP"                                      
##  [3] "GO.MF"                                      
##  [4] "GO.CC"                                      
##  [5] "JA Miller at AIBS"                          
##  [6] "Chip-X enrichment analysis (ChEA)"          
##  [7] "Brain"                                      
##  [8] "JAM"                                        
##  [9] "Prenatal brain"                             
## [10] "Brain region markers"                       
## [11] "Cortex"                                     
## [12] "Brain region marker enriched gene sets"     
## [13] "WGCNA"                                      
## [14] "BrainRegionMarkers"                         
## [15] "BrainRegionMarkers.HBA"                     
## [16] "BrainRegionMarkers.HBA.localMarker(top200)" 
## [17] "Postnatal brain"                            
## [18] "ImmunePathways"                             
## [19] "Markers of cortex layers"                   
## [20] "BrainLists"                                 
## [21] "Cell type markers"                          
## [22] "Germinal brain"                             
## [23] "BrainRegionMarkers.HBA.globalMarker(top200)"
## [24] "Accelerated evolution"                      
## [25] "Postmitotic brain"                          
## [26] "BrainLists.Blalock_AD"                      
## [27] "BrainLists.DiseaseGenes"                    
## [28] "BloodAtlases"                               
## [29] "Verge Disease Genes"                        
## [30] "BloodAtlases.Whitney"                       
## [31] "BrainLists.JAXdiseaseGene"                  
## [32] "BrainLists.MO"                              
## [33] "Age-associated genes"                       
## [34] "BrainLists.Lu_Aging"                        
## [35] "Cell type marker enriched gene sets"        
## [36] "BrainLists.CA1vsCA3"                        
## [37] "BrainLists.MitochondrialType"               
## [38] "BrainLists.MO.2+_26Mar08"                   
## [39] "BrainLists.MO.Sugino"                       
## [40] "BloodAtlases.Gnatenko2"                     
## [41] "BloodAtlases.Kabanova"                      
## [42] "BrainLists.Voineagu"                        
## [43] "StemCellLists"                              
## [44] "StemCellLists.Lee"
# Perform Enrichment Analysis
enrichment = enrichmentAnalysis(classLabels = EA_dataset$module, identifiers = EA_dataset$entrezgene,
                                refCollection = combined_col, useBackground = 'given', 
                                threshold = 1e-4, thresholdType = 'Bonferroni',
                                getOverlapEntrez = FALSE, getOverlapSymbols = TRUE)
##  enrichmentAnalysis: preparing data..
##   ..working on label set 1 ..


Results


kable(enrichment$enrichmentTable %>% filter(class==top_modules[1]) %>% 
      dplyr::select(dataSetID, shortDataSetName, inGroups, Bonferroni, FDR,
                    effectiveClassSize, effectiveSetSize, nCommonGenes),
      caption = paste0('Enriched terms for module ', top_modules[1], ' (MTcor = ',
                       round(genes_info$MTcor[genes_info$Module==top_modules[1]][1],4), ')'))
Enriched terms for module #F464E5 (MTcor = 0.9999)
dataSetID shortDataSetName inGroups Bonferroni FDR effectiveClassSize effectiveSetSize nCommonGenes
JAMiller.AIBS.000434 Genes bound by RELA in HUMAN FIBROSARCOMA from PMID 24523406 JA Miller at AIBS|Chip-X enrichment analysis (ChEA) 0.0000002 0.0e+00 642 1043 89
GO:0001525 angiogenesis GO|GO.BP 0.0000071 1.0e-07 642 434 48
GO:0072358 cardiovascular system development GO|GO.BP 0.0000273 3.0e-07 642 624 59
GO:0007166 cell surface receptor signaling pathway GO|GO.BP 0.0000365 3.0e-07 642 2498 157
GO:0001944 vasculature development GO|GO.BP 0.0000429 4.0e-07 642 615 58
GO:0072359 circulatory system development GO|GO.BP 0.0000488 4.0e-07 642 941 77
GO:0001568 blood vessel development GO|GO.BP 0.0001549 1.3e-06 642 587 55
GO:0048514 blood vessel morphogenesis GO|GO.BP 0.0002049 1.6e-06 642 511 50
GO:0035295 tube development GO|GO.BP 0.0005272 3.7e-06 642 884 71
GO:0030198 extracellular matrix organization GO|GO.BP 0.0007586 5.1e-06 642 316 36
kable(enrichment$enrichmentTable %>% filter(class==top_modules[2]) %>% 
      dplyr::select(dataSetID, shortDataSetName, inGroups, Bonferroni, FDR,
                    effectiveClassSize, effectiveSetSize, nCommonGenes),
      caption = paste0('Enriched terms for module ', top_modules[2], ' (MTcor = ',
                       round(genes_info$MTcor[genes_info$Module==top_modules[2]][1],4), ')'))
Enriched terms for module #00BECA (MTcor = -0.9186)
dataSetID shortDataSetName inGroups Bonferroni FDR effectiveClassSize effectiveSetSize nCommonGenes
JAMiller.AIBS.000142 Highest in CP of 13-16 post-conception weeks human JA Miller at AIBS|Brain|Prenatal brain|Markers of cortex layers|Cortex|Postmitotic brain 0.00e+00 0e+00 1637 1218 301
JAMiller.AIBS.000052 CortexWGCNA 15-21 post-conception weeks C26 JA Miller at AIBS|Brain|Prenatal brain|Cortex|WGCNA 0.00e+00 0e+00 1637 726 189
JAMiller.AIBS.000150 Highest in CP of E14.5 mouse JA Miller at AIBS|Brain|Prenatal brain|Markers of cortex layers|Cortex|Postmitotic brain 0.00e+00 0e+00 1637 1276 273
JAMiller.AIBS.000155 Lowest in VZ of E14.5 mouse JA Miller at AIBS|Brain|Prenatal brain|Markers of cortex layers|Cortex 0.00e+00 0e+00 1637 1674 317
JAMiller.AIBS.000569 WGCNA humanSpecificOlivedrab2Module frontalCtx FOXP2 JA Miller at AIBS|Brain|Postnatal brain|Cortex|WGCNA 0.00e+00 0e+00 1637 4057 599
JAM:002744 Autism_differential_expression_across_at_least_one_comparison JAM|BrainLists|BrainLists.Voineagu 0.00e+00 0e+00 1637 765 175
JAM:003016 downAD_synapticTransmission JAM|BrainLists|BrainLists.Blalock_AD 0.00e+00 0e+00 1637 89 48
JAMiller.AIBS.000141 CP enriched in E14.5 mouse JA Miller at AIBS|Brain|Prenatal brain|Markers of cortex layers|Cortex|Postmitotic brain 0.00e+00 0e+00 1637 569 139
JAMiller.AIBS.000506 Genes bound by SUZ12 in MOUSE MESC from PMID 20075857 JA Miller at AIBS|Chip-X enrichment analysis (ChEA) 0.00e+00 0e+00 1637 3380 497
GO:0097458 neuron part GO|GO.CC 0.00e+00 0e+00 1637 1607 281
JAMiller.AIBS.000570 WGCNA Olivedrab2ModuleGenes with enriched ELAVL2 targets JA Miller at AIBS|Brain|Postnatal brain|Cortex|WGCNA 0.00e+00 0e+00 1637 483 119
GO:0045202 synapse GO|GO.CC 0.00e+00 0e+00 1637 1114 208
GO:0099536 synaptic signaling GO|GO.BP 0.00e+00 0e+00 1637 652 142
GO:0044456 synapse part GO|GO.CC 0.00e+00 0e+00 1637 893 175
GO:0007268 chemical synaptic transmission GO|GO.BP 0.00e+00 0e+00 1637 638 138
GO:0098916 anterograde trans-synaptic signaling GO|GO.BP 0.00e+00 0e+00 1637 638 138
GO:0099537 trans-synaptic signaling GO|GO.BP 0.00e+00 0e+00 1637 646 139
JAMiller.AIBS.000005 CPi markers at 21 post-conception weeks JA Miller at AIBS|Brain|Prenatal brain|Cortex|Markers of cortex layers|Postmitotic brain 0.00e+00 0e+00 1637 312 85
JAM:003072 Tail of Caudate Nucleus_IN_Striatum JAM|BrainRegionMarkers|BrainRegionMarkers.HBA|BrainRegionMarkers.HBA.localMarker(top200)|Brain region markers|Brain region marker enriched gene sets 0.00e+00 0e+00 1637 162 57
JAM:002805 Cerebral Cortex JAM|BrainRegionMarkers|BrainRegionMarkers.HBA|BrainRegionMarkers.HBA.globalMarker(top200)|Brain region markers|Brain region marker enriched gene sets 0.00e+00 0e+00 1637 161 56
JAMiller.AIBS.000149 Lowest in VZ of 13-16 post-conception weeks human JA Miller at AIBS|Brain|Prenatal brain|Markers of cortex layers|Cortex 0.00e+00 0e+00 1637 618 131
JAMiller.AIBS.000123 HippocampusWGCNA turquoise DGenriched upAge JA Miller at AIBS|Brain|Postnatal brain|WGCNA 0.00e+00 0e+00 1637 1104 198
JAM:002751 Basal Pons JAM|BrainRegionMarkers|BrainRegionMarkers.HBA|BrainRegionMarkers.HBA.globalMarker(top200)|Brain region markers|Brain region marker enriched gene sets 0.00e+00 0e+00 1637 167 56
GO:0043005 neuron projection GO|GO.CC 0.00e+00 0e+00 1637 1214 210
JAM:003054 subiculum_IN_Hippocampal Formation JAM|BrainRegionMarkers|BrainRegionMarkers.HBA|BrainRegionMarkers.HBA.localMarker(top200)|Brain region markers|Brain region marker enriched gene sets 0.00e+00 0e+00 1637 164 53
JAM:002739 arcuate nucleus of medulla_IN_Myelencephalon JAM|BrainRegionMarkers|BrainRegionMarkers.HBA|BrainRegionMarkers.HBA.localMarker(top200)|Brain region markers|Brain region marker enriched gene sets 0.00e+00 0e+00 1637 168 53
GO:0097060 synaptic membrane GO|GO.CC 0.00e+00 0e+00 1637 412 94
JAM:002907 inferior olivary complex_IN_Myelencephalon JAM|BrainRegionMarkers|BrainRegionMarkers.HBA|BrainRegionMarkers.HBA.localMarker(top200)|Brain region markers|Brain region marker enriched gene sets 0.00e+00 0e+00 1637 170 53
JAMiller.AIBS.000504 Genes bound by SUZ12 in mouse MESC from PMID 18692474 JA Miller at AIBS|Chip-X enrichment analysis (ChEA) 0.00e+00 0e+00 1637 1362 222
JAM:002824 Dentate Nucleus_IN_Cerebellar Nucleus JAM|BrainRegionMarkers|BrainRegionMarkers.HBA|BrainRegionMarkers.HBA.localMarker(top200)|Brain region markers|Brain region marker enriched gene sets 0.00e+00 0e+00 1637 165 51
GO:0030424 axon GO|GO.CC 0.00e+00 0e+00 1637 574 116
GO:1902495 transmembrane transporter complex GO|GO.CC 0.00e+00 0e+00 1637 289 72
JAM:003058 Substantia Nigra, pars compacta_IN_Mesencephalon JAM|BrainRegionMarkers|BrainRegionMarkers.HBA|BrainRegionMarkers.HBA.localMarker(top200)|Brain region markers|Brain region marker enriched gene sets 0.00e+00 0e+00 1637 169 51
GO:0034702 ion channel complex GO|GO.CC 0.00e+00 0e+00 1637 269 68
JAMiller.AIBS.000364 Genes bound by MTF2 in MOUSE MESC from PMID 20144788 JA Miller at AIBS|Chip-X enrichment analysis (ChEA) 0.00e+00 0e+00 1637 2292 330
JAM:002918 lateral medullary reticular group_IN_Myelencephalon JAM|BrainRegionMarkers|BrainRegionMarkers.HBA|BrainRegionMarkers.HBA.localMarker(top200)|Brain region markers|Brain region marker enriched gene sets 0.00e+00 0e+00 1637 162 49
GO:1990351 transporter complex GO|GO.CC 1.00e-07 0e+00 1637 297 72
JAM:002764 downAging_mitochondria_synapse JAM|BrainLists|BrainLists.Lu_Aging 1.00e-07 0e+00 1637 391 86
GO:0045211 postsynaptic membrane GO|GO.CC 1.00e-07 0e+00 1637 309 73
JAMiller.AIBS.000503 Genes bound by SUZ12 in mouse MESC from PMID 18555785 JA Miller at AIBS|Chip-X enrichment analysis (ChEA) 2.00e-07 0e+00 1637 765 138
JAM:002991 downAD_synapticTransmission JAM|BrainLists|BrainLists.Blalock_AD 3.00e-07 0e+00 1637 97 35
JAM:002967 Occipital Lobe_IN_Cerebral Cortex JAM|BrainRegionMarkers|BrainRegionMarkers.HBA|BrainRegionMarkers.HBA.localMarker(top200)|Brain region markers|Brain region marker enriched gene sets 4.00e-07 0e+00 1637 149 45
GO:0050877 nervous system process GO|GO.BP 5.00e-07 0e+00 1637 886 153
JAM:002964 Nucleus Accumbens_IN_Striatum JAM|BrainRegionMarkers|BrainRegionMarkers.HBA|BrainRegionMarkers.HBA.localMarker(top200)|Brain region markers|Brain region marker enriched gene sets 6.00e-07 0e+00 1637 161 47
GO:0098794 postsynapse GO|GO.CC 6.00e-07 0e+00 1637 591 113
GO:0120025 plasma membrane bounded cell projection GO|GO.CC 8.00e-07 0e+00 1637 1919 280
GO:0034703 cation channel complex GO|GO.CC 9.00e-07 0e+00 1637 203 54
JAM:002882 Hippocampal Formation JAM|BrainRegionMarkers|BrainRegionMarkers.HBA|BrainRegionMarkers.HBA.globalMarker(top200)|Brain region markers|Brain region marker enriched gene sets 1.00e-06 0e+00 1637 169 48
JAM:002763 downAD_metalIonTransport_glycoprotein JAM|BrainLists|BrainLists.Blalock_AD 1.10e-06 0e+00 1637 283 67
JAM:003085 trochlear nucleus_IN_Mesencephalon JAM|BrainRegionMarkers|BrainRegionMarkers.HBA|BrainRegionMarkers.HBA.localMarker(top200)|Brain region markers|Brain region marker enriched gene sets 1.20e-06 0e+00 1637 153 45
GO:0007267 cell-cell signaling GO|GO.BP 1.60e-06 0e+00 1637 1439 221
GO:0042995 cell projection GO|GO.CC 2.60e-06 0e+00 1637 1989 286
GO:0050804 modulation of chemical synaptic transmission GO|GO.BP 4.40e-06 1e-07 1637 411 85
GO:0005886 plasma membrane GO|GO.CC 4.50e-06 1e-07 1637 4277 541
GO:0099177 regulation of trans-synaptic signaling GO|GO.BP 5.00e-06 1e-07 1637 412 85
GO:0098793 presynapse GO|GO.CC 6.10e-06 1e-07 1637 462 92
GO:0071944 cell periphery GO|GO.CC 7.10e-06 1e-07 1637 4381 551
GO:0098978 glutamatergic synapse GO|GO.CC 7.70e-06 1e-07 1637 341 74
JAM:002920 Lateral Nucleus_IN_Amygdala JAM|BrainRegionMarkers|BrainRegionMarkers.HBA|BrainRegionMarkers.HBA.localMarker(top200)|Brain region markers|Brain region marker enriched gene sets 8.40e-06 1e-07 1637 167 46
GO:0005216 ion channel activity GO|GO.MF 8.80e-06 1e-07 1637 362 77
JAM:002769 downAD_mitochondrion JAM|BrainLists|BrainLists.Blalock_AD 1.04e-05 1e-07 1637 265 62
JAMiller.AIBS.000463 Genes bound by SMAD4 in HUMAN A2780 from PMID 21799915 JA Miller at AIBS|Chip-X enrichment analysis (ChEA) 1.18e-05 1e-07 1637 2085 294
JAMiller.AIBS.000042 CortexWGCNA 15-21 post-conception weeks C16 SPenriched JA Miller at AIBS|Brain|Prenatal brain|Cortex|WGCNA 1.21e-05 1e-07 1637 198 51
GO:0036477 somatodendritic compartment GO|GO.CC 1.63e-05 2e-07 1637 795 136
GO:0044463 cell projection part GO|GO.CC 1.91e-05 2e-07 1637 1345 205
GO:0120038 plasma membrane bounded cell projection part GO|GO.CC 1.91e-05 2e-07 1637 1345 205
GO:0006812 cation transport GO|GO.BP 1.93e-05 2e-07 1637 1007 163
GO:0005261 cation channel activity GO|GO.MF 2.04e-05 2e-07 1637 282 64
GO:0022803 passive transmembrane transporter activity GO|GO.MF 2.14e-05 2e-07 1637 389 80
GO:0022838 substrate-specific channel activity GO|GO.MF 2.96e-05 3e-07 1637 371 77
GO:0003008 system process GO|GO.BP 3.92e-05 4e-07 1637 1505 223
JAM:002866 GlutamatergicNeuronsInMouseCortex_Sugino JAM|BrainLists|BrainLists.MO|BrainLists.MO.Sugino|Cell type marker enriched gene sets|Brain region marker enriched gene sets 4.29e-05 4e-07 1637 367 76
GO:0015267 channel activity GO|GO.MF 4.53e-05 4e-07 1637 388 79
JAMiller.AIBS.000436 Genes bound by REST in MOUSE MESC from PMID 21632747 JA Miller at AIBS|Chip-X enrichment analysis (ChEA) 4.54e-05 4e-07 1637 1675 243
JAMiller.AIBS.000505 Genes bound by SUZ12 in MOUSE MESC from PMID 18974828 JA Miller at AIBS|Chip-X enrichment analysis (ChEA) 5.82e-05 5e-07 1637 1387 208
GO:0043269 regulation of ion transport GO|GO.BP 5.94e-05 5e-07 1637 613 110

Save Enrichment Analysis results

save(enrichment, file='./../Data/enrichmentAnalysis.RData')



Session info

sessionInfo()
## R version 3.6.0 (2019-04-26)
## Platform: x86_64-redhat-linux-gnu (64-bit)
## Running under: Scientific Linux 7.6 (Nitrogen)
## 
## Matrix products: default
## BLAS/LAPACK: /usr/lib64/R/lib/libRblas.so
## 
## locale:
##  [1] LC_CTYPE=en_GB.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_GB.UTF-8        LC_COLLATE=en_GB.UTF-8    
##  [5] LC_MONETARY=en_GB.UTF-8    LC_MESSAGES=en_GB.UTF-8   
##  [7] LC_PAPER=en_GB.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_GB.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats4    parallel  stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] org.Hs.eg.db_3.10.0                      
##  [2] BrainDiseaseCollection_1.00              
##  [3] anRichment_1.01-2                        
##  [4] TxDb.Mmusculus.UCSC.mm10.knownGene_3.10.0
##  [5] TxDb.Hsapiens.UCSC.hg19.knownGene_3.2.2  
##  [6] GenomicFeatures_1.38.2                   
##  [7] GenomicRanges_1.38.0                     
##  [8] GenomeInfoDb_1.22.0                      
##  [9] anRichmentMethods_0.90-1                 
## [10] WGCNA_1.68                               
## [11] fastcluster_1.1.25                       
## [12] dynamicTreeCut_1.63-1                    
## [13] GO.db_3.10.0                             
## [14] AnnotationDbi_1.48.0                     
## [15] IRanges_2.20.2                           
## [16] S4Vectors_0.24.3                         
## [17] Biobase_2.46.0                           
## [18] BiocGenerics_0.32.0                      
## [19] biomaRt_2.42.0                           
## [20] knitr_1.24                               
## [21] doParallel_1.0.15                        
## [22] iterators_1.0.12                         
## [23] foreach_1.4.7                            
## [24] polycor_0.7-10                           
## [25] expss_0.10.1                             
## [26] GGally_1.4.0                             
## [27] gridExtra_2.3                            
## [28] viridis_0.5.1                            
## [29] viridisLite_0.3.0                        
## [30] RColorBrewer_1.1-2                       
## [31] dendextend_1.13.3                        
## [32] plotly_4.9.2                             
## [33] glue_1.3.1                               
## [34] reshape2_1.4.3                           
## [35] forcats_0.4.0                            
## [36] stringr_1.4.0                            
## [37] dplyr_0.8.3                              
## [38] purrr_0.3.3                              
## [39] readr_1.3.1                              
## [40] tidyr_1.0.2                              
## [41] tibble_2.1.3                             
## [42] ggplot2_3.2.1                            
## [43] tidyverse_1.3.0                          
## 
## loaded via a namespace (and not attached):
##   [1] readxl_1.3.1                backports_1.1.5            
##   [3] Hmisc_4.2-0                 BiocFileCache_1.10.2       
##   [5] plyr_1.8.5                  lazyeval_0.2.2             
##   [7] splines_3.6.0               crosstalk_1.0.0            
##   [9] BiocParallel_1.20.1         robust_0.4-18.2            
##  [11] digest_0.6.24               htmltools_0.4.0            
##  [13] fansi_0.4.1                 magrittr_1.5               
##  [15] checkmate_1.9.4             memoise_1.1.0              
##  [17] fit.models_0.5-14           cluster_2.0.8              
##  [19] annotate_1.64.0             Biostrings_2.54.0          
##  [21] modelr_0.1.5                matrixStats_0.55.0         
##  [23] askpass_1.1                 prettyunits_1.0.2          
##  [25] colorspace_1.4-1            blob_1.2.1                 
##  [27] rvest_0.3.5                 rappdirs_0.3.1             
##  [29] rrcov_1.4-7                 haven_2.2.0                
##  [31] xfun_0.8                    crayon_1.3.4               
##  [33] RCurl_1.95-4.12             jsonlite_1.6               
##  [35] genefilter_1.68.0           impute_1.60.0              
##  [37] survival_2.44-1.1           gtable_0.3.0               
##  [39] zlibbioc_1.32.0             XVector_0.26.0             
##  [41] DelayedArray_0.12.2         DEoptimR_1.0-8             
##  [43] scales_1.1.0                mvtnorm_1.0-11             
##  [45] DBI_1.1.0                   Rcpp_1.0.3                 
##  [47] xtable_1.8-4                progress_1.2.2             
##  [49] htmlTable_1.13.1            foreign_0.8-71             
##  [51] bit_1.1-15.2                preprocessCore_1.48.0      
##  [53] Formula_1.2-3               htmlwidgets_1.5.1          
##  [55] httr_1.4.1                  ellipsis_0.3.0             
##  [57] acepack_1.4.1               farver_2.0.3               
##  [59] pkgconfig_2.0.3             reshape_0.8.8              
##  [61] XML_3.99-0.3                nnet_7.3-12                
##  [63] dbplyr_1.4.2                locfit_1.5-9.1             
##  [65] later_1.0.0                 labeling_0.3               
##  [67] tidyselect_0.2.5            rlang_0.4.4                
##  [69] munsell_0.5.0               cellranger_1.1.0           
##  [71] tools_3.6.0                 cli_2.0.1                  
##  [73] generics_0.0.2              RSQLite_2.2.0              
##  [75] broom_0.5.4                 fastmap_1.0.1              
##  [77] evaluate_0.14               yaml_2.2.0                 
##  [79] bit64_0.9-7                 fs_1.3.1                   
##  [81] robustbase_0.93-5           nlme_3.1-139               
##  [83] mime_0.9                    xml2_1.2.2                 
##  [85] compiler_3.6.0              rstudioapi_0.10            
##  [87] curl_4.3                    reprex_0.3.0               
##  [89] geneplotter_1.64.0          pcaPP_1.9-73               
##  [91] stringi_1.4.6               highr_0.8                  
##  [93] lattice_0.20-38             Matrix_1.2-17              
##  [95] vctrs_0.2.2                 pillar_1.4.3               
##  [97] lifecycle_0.1.0             data.table_1.12.8          
##  [99] bitops_1.0-6                httpuv_1.5.2               
## [101] rtracklayer_1.46.0          R6_2.4.1                   
## [103] latticeExtra_0.6-28         promises_1.1.0             
## [105] codetools_0.2-16            MASS_7.3-51.4              
## [107] assertthat_0.2.1            SummarizedExperiment_1.16.1
## [109] DESeq2_1.26.0               openssl_1.4.1              
## [111] withr_2.1.2                 GenomicAlignments_1.22.1   
## [113] Rsamtools_2.2.2             GenomeInfoDbData_1.2.2     
## [115] hms_0.5.3                   grid_3.6.0                 
## [117] rpart_4.1-15                rmarkdown_1.14             
## [119] Cairo_1.5-10                shiny_1.4.0                
## [121] lubridate_1.7.4             base64enc_0.1-3